Context

Jaeger chez Jessberger provided a Seurat object with clustering results/genotypes/batch effects and asked for a differential expression analysis with batch as blocking factors (Friday 16th Aug 2019).

To be more precise, Ascl1_5day_2 and Ascl1_12wk_2 show batch effects they would like to subtract. Analysis are run using edgeR’s glmQLFit with batch as blocking factor; and logistic regression (LR) with batch as latent variables.

Ascl vs Gli comparisons are relevant biologically; days (5d, 12 week etc) are unimportant.

Only genes with at least TPM above 1 in more than 25% of the cells are included into the edgeR glmQLFit to avoid too many false positives for lowly expressed genes.

Only genes represented in 20% of the cells of at least one group are included into the logistic regression analysis.

Datatables can be downloaded to be browsed using R or spreadsheet software.

Some of the markers are plotted as ridgeplots, but the comprehensive list of significant results (adj pvalue < 0.05) can be browsed/downloaded as datatables.

P-value adjustment Bonferroni for LR tests, FDR for edgeR’s glmQLFit.

Last results are cluster-wise, e.g. cells are compared according to both genotype (Ascl or Gli) and cluster (MN, IN, etc; annotated by Jessberger people) they belong to. We provide two analysis for that:

-First markers that are specific for each one of the genotype:cluster combination (that allows to distinguish them from the others). Instead of comparing Ascl vs Gli1 in cluster 1, then in cluster 2 etc what I did was to get the markers that distinguish Ascl:cluster1 from everything else, Ascl:cluster2 from everything else, etc. Thus keeping the ‘background’ in all comparisons. - Second, genotypes comparison within clusters (only mn_ascl vs mn_gli, in_ascl vs in_gli etc).

Load

analysis_tag <- '03_diff_expression'
library(pheatmap)
library(knitr)
library(scater)
library(Seurat)
library(Cairo)
library(edgeR)
library(data.table)
## library(biomaRt)
library(dplyr)
## library(reshape2)
#library(iSEE)
## library(Matrix)
## library(scran)
## library(mvoutlier)
## library(shiny)
## library(kableExtra)
## library(velocyto.R)
## library(cowplot)
## library(corrplot)
library(DT)
library(devtools)
## library(viridis)
## library("CellMixS")
library(gplots)
library(pheatmap)
ac <- function(col, alpha = 1) {
    
    apply(sapply(col, col2rgb)/255, 2, function(x) rgb(x[1], 
        x[2], x[3], alpha = alpha))
}

##' subsets a dge to cells annotated as `genotype` in column
## ' `clusterwise`; 25% cells with over 1 TPM to be returned
##'
##' 
##' @title dgelist subset
##' @param dge dgelist
##' @param genotypes charvector
##' @param tpm_mat tpm matrix
##' @return a dgelist object
##' @author izaskun.mallona@gmail.com
subset_cells <- function(dge, genotypes, tpm_mat) {
    mod <- dge
    cells <- rownames(dge$samples[dge$samples$clusterwise %in% 
        genotypes, ])
    
    keep <- rownames(tpm_mat)[rowSums(tpm_mat[, cells] > 1) >= 
        0.25 * ncol(dge$counts[, cells])]
    
    mod$counts <- dge$counts[keep, cells]
    mod$samples <- dge$samples[cells, ]
    return(mod)
}

Data load

Provided by Jaeger

d <- readRDS(file.path("..", "data", "Seurat_object_Gli_Ascl_combi_meta.Robj"))

Checking metadata

table(d@meta.data$group, d@meta.data$cluster)
            
             ndNSC dNSC  IN  MN
  Ascl1_5day    35  109 172 185
  Ascl1_12wk    30   26  65 112
  Gli1_5day     63   25  19  28
  Gli1_12wk     20   19  37 110
table(d@meta.data$group, d@meta.data$batch_group)
            
             Ascl1_12wk_1 Ascl1_12wk_2 Ascl1_5day_1 Ascl1_5day_2 Gli1_12wk
  Ascl1_5day            0            0          264          237         0
  Ascl1_12wk           57          176            0            0         0
  Gli1_5day             0            0            0            0         0
  Gli1_12wk             0            0            0            0       186
            
             Gli1_5day_1 Gli1_5day_2
  Ascl1_5day           0           0
  Ascl1_12wk           0           0
  Gli1_5day           55          80
  Gli1_12wk            0           0

Diff expression All cells (all clusters) Gli1 vs Ascl1

edgeR DGEList object (full data)

dge_full <- DGEList(counts = as.matrix(GetAssayData(d@assays$RNA, 
    slot = "counts")), samples = d@meta.data)

Annotating genes (Asia biomart mirror because european was down while rendering)

tmp <- list()

tmp$ids = rownames(dge_full$counts)

tmp$filters = "mgi_symbol"
tmp$attributes = c(tmp$filters, "ensembl_gene_id", "chromosome_name", 
    "gene_biotype", "start_position", "end_position")
tmp$biomart = "ENSEMBL_MART_ENSEMBL"
tmp$dataset = "mmusculus_gene_ensembl"
tmp$host = "asia.ensembl.org"
tmp$bmart <- biomaRt::useMart(biomart = tmp$biomart, dataset = tmp$dataset, 
    host = tmp$host)

feature_info <- biomaRt::getBM(attributes = tmp$attributes, filters = tmp$filters, 
    values = tmp$ids, mart = tmp$bmart)
Batch submitting query [=>------------------------] 6% eta: 18s
Batch submitting query [=>------------------------] 9% eta: 23s
Batch submitting query [==>-----------------------] 12% eta: 25s
Batch submitting query [===>----------------------] 16% eta: 26s
Batch submitting query [====>---------------------] 19% eta: 26s
Batch submitting query [=====>--------------------] 22% eta: 26s
Batch submitting query [=====>--------------------] 25% eta: 25s
Batch submitting query [======>-------------------] 28% eta: 24s
Batch submitting query [=======>------------------] 31% eta: 24s
Batch submitting query [========>-----------------] 34% eta: 23s
Batch submitting query [=========>----------------] 38% eta: 22s
Batch submitting query [==========>---------------] 41% eta: 21s
Batch submitting query [==========>---------------] 44% eta: 20s
Batch submitting query [===========>--------------] 47% eta: 19s
Batch submitting query [============>-------------] 50% eta: 18s
Batch submitting query [=============>------------] 53% eta: 17s
Batch submitting query [==============>-----------] 56% eta: 16s
Batch submitting query [==============>-----------] 59% eta: 15s
Batch submitting query [===============>----------] 62% eta: 14s
Batch submitting query [================>---------] 66% eta: 13s
Batch submitting query [=================>--------] 69% eta: 12s
Batch submitting query [==================>-------] 72% eta: 10s
Batch submitting query [===================>------] 75% eta: 9s
Batch submitting query [===================>------] 78% eta: 8s
Batch submitting query [====================>-----] 81% eta: 7s
Batch submitting query [=====================>----] 84% eta: 6s
Batch submitting query [======================>---] 88% eta: 5s
Batch submitting query [=======================>--] 91% eta: 4s Batch
submitting query [=======================>--] 94% eta: 2s Batch submitting
query [========================>-] 97% eta: 1s Batch submitting query
[==========================] 100% eta: 0s
mm <- match(tmp$ids, feature_info[[tmp$filters]])
feature_info_full <- feature_info[mm, ]
old_rdata <- data.frame(id = rownames(dge_full$counts))
keep_cols <- !(colnames(old_rdata) %in% colnames(feature_info_full))
gene_meta_data <- cbind(old_rdata[, keep_cols], feature_info_full)
rownames(gene_meta_data) <- gene_meta_data[, 1]
gene_meta_data$length <- gene_meta_data$end_position - gene_meta_data$start_position
## rm(tmp, old_rdata, keep_cols, mm, feature_info)

Removing mitochondrial genes

keep <- !grepl("mt-", rownames(dge_full))
table(keep)
keep
FALSE  TRUE 
   13 15502 
dge_full <- dge_full[keep, ]
DefaultAssay(d) <- "RNA"
keep <- !grepl("mt-", rownames(d@assays$RNA))
table(keep)
keep
FALSE  TRUE 
   13 15502 
d <- d[keep, ]

Generating unexact TPMs from raw counts (because we take gene lengths and didn’t quantify isoforms)

## removing NAs

gene_meta_data <- na.omit(gene_meta_data)

genes <- intersect(rownames(dge_full$counts), rownames(gene_meta_data))

dge_full$counts <- dge_full$counts[genes, ]

x <- dge_full$counts/gene_meta_data[rownames(dge_full$counts), 
    "length"]
tpm_mat <- t(t(x) * 1e+06/colSums(x))

Filtering in only genes with an estimated TPM above 1 in more than 25% of the cells, as in https://www.nature.com/articles/nmeth.4612. TPMs because of the sequencing length bias https://f1000research.com/articles/6-595.

dge <- dge_full
keep <- rowSums(tpm_mat > 1) >= 0.25 * ncol(dge$counts)
table(keep)
keep
FALSE  TRUE 
10235  4895 
dge <- dge[keep, ]
batch <- dge$samples$batch_group
genotype <- dge$samples$type
cluster <- dge$samples$cluster

Batch and cluster are full rank, trying to split the batch into a meaningful design

alt_batch <- as.character(batch)

alt_batch[alt_batch == "Ascl1_12wk_1"] <- "B1"
alt_batch[alt_batch == "Ascl1_12wk_2"] <- "B2"
alt_batch[alt_batch == "Ascl1_5day_1"] <- "B1"
alt_batch[alt_batch == "Ascl1_5day_2"] <- "B2"
alt_batch[alt_batch == "Gli1_12wk"] <- "B1"
alt_batch[alt_batch == "Gli1_5day_1"] <- "B1"
alt_batch[alt_batch == "Gli1_5day_2"] <- "B1"
alt_batch <- as.factor(alt_batch)
table(alt_batch)
alt_batch
 B1  B2 
642 413 
design <- model.matrix(~0 + genotype + alt_batch)
# norm factors
dge <- calcNormFactors(dge)

# Fit the NB GLMs with QL methods
dge <- estimateDisp(dge, design)
fit <- glmQLFit(dge, design)
limma::plotMDS(dge, col = as.numeric(as.factor(dge$samples$type)), 
    pch = 19)

plotBCV(dge)

plotQLDisp(fit)

Comparing Gli and Ascl genotypes with glmQLF

avsgli <- makeContrasts(genotypeAscl1 - genotypeGli1, levels = design)
qlf_genot <- glmQLFTest(fit, contrast = avsgli)
tt_genot <- topTags(qlf_genot, adjust.method = "BH", sort.by = "PValue", 
    n = 1e+06, p.value = 0.05)
tt_genot_df <- data.frame(gene = rownames(tt_genot), tt_genot)
tt_genot_df$analysis <- "ascl_vs_gli_glmqlf"

DT::datatable(tt_genot_df %>% as.data.frame() %>% dplyr::mutate_if(is.numeric, 
    funs(round(., 2))), extensions = c("Buttons", "FixedColumns"), 
    rownames = FALSE, options = list(dom = "Bfrtip", scrollX = TRUE, 
        fixedColumns = list(leftColumns = 1), buttons = c("csv", 
            "excel")))
Warning: funs() is soft deprecated as of dplyr 0.8.0
Please use a list of either functions or lambdas: 

  # Simple named list: 
  list(mean = mean, median = median)

  # Auto named with `tibble::lst()`: 
  tibble::lst(mean, median)

  # Using lambdas
  list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
This warning is displayed once per session.
hist(tt_genot$table$PValue, 50)

hist(tt_genot$table$FDR, 50)

plotSmear(qlf_genot)

Heat map of top 50 differentially expressed genes, logCPM as in https://f1000research.com/articles/5-1438

logCPM <- cpm(dge, prior.count = 2, log = TRUE)
rownames(logCPM) <- rownames(dge$counts)
colnames(logCPM) <- rownames(dge$samples)
## table is already sorted
o <- head(rownames(tt_genot$table), 50)

logCPM <- logCPM[o, ]
logCPM <- t(scale(t(logCPM)))

col.pan <- colorpanel(100, "blue", "white", "red")

## sampling only 100 cells
set.seed(1)
logCPM <- logCPM[, sample(1:ncol(logCPM), 100)]



pheatmap(logCPM, color = col.pan, annotation_col = dge$samples[, 
    c("cluster", "batch_group", "type")])
DefaultAssay(d) <- "RNA"
Idents(d) <- d@meta.data$combi
RidgePlot(d, head(tt_genot_df$gene), ncol = 2)

Comparing Gli and Ascl genotypes with logistic regression

In this case, we require 20% of the cells in either Ascl1 or Gli1 to express the gene.

Note about differential expression results: it hasthe following columns, - gene: the gene symbol - p_val : p_val (unadjusted) - avg_logFC : log fold-chage of the average expression between the two groups. Positive values indicate that the feature is more highly expressed in the first group. - pct.1 : The percentage of cells where the feature is detected in the first group - pct.2 : The percentage of cells where the feature is detected in the second group - p_val_adj : Adjusted p-value, based on bonferroni correction using all features in the dataset. - analysis: a tag describing the analysis carried out

alt_batch <- as.character(d@meta.data$batch)


alt_batch[alt_batch == "Ascl1_12wk_1"] <- "B1"
alt_batch[alt_batch == "Ascl1_12wk_2"] <- "B2"
alt_batch[alt_batch == "Ascl1_5day_1"] <- "B1"
alt_batch[alt_batch == "Ascl1_5day_2"] <- "B2"
alt_batch[alt_batch == "Gli1_12wk"] <- "B1"
alt_batch[alt_batch == "Gli1_5day_1"] <- "B1"
alt_batch[alt_batch == "Gli1_5day_2"] <- "B1"
alt_batch <- as.factor(alt_batch)
table(alt_batch)
alt_batch
 B1  B2 
642 413 
d@meta.data$alt_batch <- alt_batch


DefaultAssay(d) <- "RNA"
Idents(d) <- "type"

markers <- FindMarkers(d, ident.1 = "Ascl1", ident.2 = "Gli1", 
    test.use = "LR", latent.vars = "alt_batch", min.pct = 0.2)
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
markers <- markers[order(markers$p_val_adj), ]
markers <- markers[markers$p_val_adj < 0.05, ]
markers$analysis <- "ascl_vs_gli_lr"
markers <- data.frame(gene = rownames(markers), markers)

DT::datatable(markers %>% as.data.frame() %>% dplyr::mutate_if(is.numeric, 
    funs(round(., 2))), extensions = c("Buttons", "FixedColumns"), 
    rownames = FALSE, options = list(dom = "Bfrtip", scrollX = TRUE, 
        fixedColumns = list(leftColumns = 1), buttons = c("csv", 
            "excel")))
Idents(d) <- d@meta.data$combi
RidgePlot(d, head(rownames(markers)), ncol = 2)

Comparing the overlap between LR and edgeR

edge <- data.frame(gene = rownames(tt_genot$table), tt_genot$table)

lr <- data.frame(gene = rownames(markers), markers)

merged <- merge(edge, lr, by = "gene", all = TRUE)

nrow(edge)
[1] 1335
nrow(lr)
[1] 320
length(intersect(edge$gene, lr$gene))
[1] 264
merged$logFDR <- log10(merged$FDR)
merged$logp_val_adj <- log10(merged$p_val_adj)

plot(merged[c("logFC", "avg_logFC", "logFDR", "logp_val_adj")])

Unsupervised genotype:cluster-wise DE analysis

Getting the markers that are characteristic of each genotype:cluster combination (as compared to the rest).

DefaultAssay(d) <- "RNA"
Idents(d) <- "combi"
markers_combi <- FindAllMarkers(d, test.use = "LR", latent.vars = "alt_batch", 
    min.pct = 0.2)
Calculating cluster Gli1_ndNSC
Calculating cluster Ascl1_ndNSC
Calculating cluster Gli1_dNSC
Calculating cluster Ascl1_dNSC
Calculating cluster Gli1_IN
Calculating cluster Ascl1_IN
Calculating cluster Gli1_MN
Calculating cluster Ascl1_MN
markers_combi <- markers_combi[markers_combi$p_val_adj < 0.05, 
    ]
markers_combi$analysis <- "all_markers_lr"

DT::datatable(markers_combi %>% as.data.frame() %>% dplyr::mutate_if(is.numeric, 
    funs(round(., 2))), extensions = c("Buttons", "FixedColumns"), 
    rownames = FALSE, options = list(dom = "Bfrtip", scrollX = TRUE, 
        fixedColumns = list(leftColumns = 1), buttons = c("csv", 
            "excel")))

Planned genotype:cluster-wise comparisons

Finally, genotypes comparison within clusters, e.g. for cluster MN comparing cells from mn_ascl vs mn_gli; for IN cells, coparing in_ascl vs in_gli etc.

Ascl vs Gli for cluster ndNSC

DefaultAssay(d) <- "RNA"
Idents(d) <- "combi"

markers <- FindMarkers(d, ident.1 = "Ascl1_ndNSC", ident.2 = "Gli1_ndNSC", 
    test.use = "LR", latent.vars = "alt_batch", min.pct = 0.2)

markers <- markers[order(markers$p_val_adj), ]
markers <- markers[markers$p_val_adj < 0.05, ]
markers$analysis <- "asclndnsc_glindnsc_lr"
markers <- data.frame(gene = rownames(markers), markers)

DT::datatable(markers %>% as.data.frame() %>% dplyr::mutate_if(is.numeric, 
    funs(round(., 2))), extensions = c("Buttons", "FixedColumns"), 
    rownames = FALSE, options = list(dom = "Bfrtip", scrollX = TRUE, 
        fixedColumns = list(leftColumns = 1), buttons = c("csv", 
            "excel")))
Idents(d) <- "combi"
RidgePlot(d, head(rownames(markers)), ncol = 2)

Ascl vs Gli for cluster dNSC

DefaultAssay(d) <- "RNA"
Idents(d) <- "combi"

markers <- FindMarkers(d, ident.1 = "Ascl1_dNSC", ident.2 = "Gli1_dNSC", 
    test.use = "LR", latent.vars = "alt_batch", min.pct = 0.2)

markers <- markers[order(markers$p_val_adj), ]
markers <- markers[markers$p_val_adj < 0.05, ]
markers$analysis <- "ascldnsc_glidnsc_lr"
markers <- data.frame(gene = rownames(markers), markers)

DT::datatable(markers %>% as.data.frame() %>% dplyr::mutate_if(is.numeric, 
    funs(round(., 2))), extensions = c("Buttons", "FixedColumns"), 
    rownames = FALSE, options = list(dom = "Bfrtip", scrollX = TRUE, 
        fixedColumns = list(leftColumns = 1), buttons = c("csv", 
            "excel")))
Idents(d) <- "combi"
RidgePlot(d, head(rownames(markers)), ncol = 2)

Ascl vs Gli for cluster IN

DefaultAssay(d) <- "RNA"
Idents(d) <- "combi"

markers <- FindMarkers(d, ident.1 = "Ascl1_IN", ident.2 = "Gli1_IN", 
    test.use = "LR", latent.vars = "alt_batch", min.pct = 0.2)

markers <- markers[order(markers$p_val_adj), ]
markers <- markers[markers$p_val_adj < 0.05, ]
markers$analysis <- "asclin_gliin_lr"
markers <- data.frame(gene = rownames(markers), markers)

DT::datatable(markers %>% as.data.frame() %>% dplyr::mutate_if(is.numeric, 
    funs(round(., 2))), extensions = c("Buttons", "FixedColumns"), 
    rownames = FALSE, options = list(dom = "Bfrtip", scrollX = TRUE, 
        fixedColumns = list(leftColumns = 1), buttons = c("csv", 
            "excel")))
Idents(d) <- "combi"
RidgePlot(d, head(rownames(markers)), ncol = 2)

Ascl vs Gli for cluster MN

DefaultAssay(d) <- "RNA"
Idents(d) <- "combi"

markers <- FindMarkers(d, ident.1 = "Ascl1_MN", ident.2 = "Gli1_MN", 
    test.use = "LR", latent.vars = "alt_batch", min.pct = 0.2)

markers <- markers[order(markers$p_val_adj), ]
markers <- markers[markers$p_val_adj < 0.05, ]
markers$analysis <- "asclmn_glimn_lr"
markers <- data.frame(gene = rownames(markers), markers)

DT::datatable(markers %>% as.data.frame() %>% dplyr::mutate_if(is.numeric, 
    funs(round(., 2))), extensions = c("Buttons", "FixedColumns"), 
    rownames = FALSE, options = list(dom = "Bfrtip", scrollX = TRUE, 
        fixedColumns = list(leftColumns = 1), buttons = c("csv", 
            "excel")))
Idents(d) <- "combi"
RidgePlot(d, head(rownames(markers)), ncol = 2)

Session

date()
[1] "Mon Aug 19 16:29:50 2019"
devtools::session_info()
─ Session info ──────────────────────────────────────────────────────────
 setting  value                       
 version  R version 3.6.1 (2019-07-05)
 os       Ubuntu 16.04.6 LTS          
 system   x86_64, linux-gnu           
 ui       X11                         
 language en_US                       
 collate  en_US.UTF-8                 
 ctype    en_US.UTF-8                 
 tz       Europe/Madrid               
 date     2019-08-19                  

─ Packages ──────────────────────────────────────────────────────────────
 package              * version   date       lib source        
 AnnotationDbi          1.46.0    2019-05-02 [1] Bioconductor  
 ape                    5.3       2019-03-17 [1] CRAN (R 3.6.1)
 assertthat             0.2.1     2019-03-21 [1] CRAN (R 3.6.1)
 backports              1.1.4     2019-04-10 [1] CRAN (R 3.6.1)
 beeswarm               0.2.3     2016-04-25 [1] CRAN (R 3.6.1)
 bibtex                 0.4.2     2017-06-30 [1] CRAN (R 3.6.1)
 Biobase              * 2.44.0    2019-05-02 [1] Bioconductor  
 BiocGenerics         * 0.30.0    2019-05-02 [1] Bioconductor  
 BiocNeighbors          1.2.0     2019-05-02 [1] Bioconductor  
 BiocParallel         * 1.18.0    2019-05-03 [1] Bioconductor  
 BiocSingular           1.0.0     2019-05-02 [1] Bioconductor  
 biomaRt                2.40.1    2019-06-27 [1] Bioconductor  
 bit                    1.1-14    2018-05-29 [1] CRAN (R 3.6.1)
 bit64                  0.9-7     2017-05-08 [1] CRAN (R 3.6.1)
 bitops                 1.0-6     2013-08-17 [1] CRAN (R 3.6.1)
 blob                   1.2.0     2019-07-09 [1] CRAN (R 3.6.1)
 Cairo                * 1.5-10    2019-03-28 [1] CRAN (R 3.6.1)
 callr                  3.3.1     2019-07-18 [1] CRAN (R 3.6.1)
 caTools                1.17.1.2  2019-03-06 [1] CRAN (R 3.6.1)
 cli                    1.1.0     2019-03-19 [1] CRAN (R 3.6.1)
 cluster                2.1.0     2019-06-19 [1] CRAN (R 3.6.1)
 codetools              0.2-16    2018-12-24 [1] CRAN (R 3.6.1)
 colorspace             1.4-1     2019-03-18 [1] CRAN (R 3.6.1)
 cowplot                1.0.0     2019-07-11 [1] CRAN (R 3.6.1)
 crayon                 1.3.4     2017-09-16 [1] CRAN (R 3.6.1)
 crosstalk              1.0.0     2016-12-21 [1] CRAN (R 3.6.1)
 curl                   3.3       2019-01-10 [1] CRAN (R 3.6.1)
 data.table           * 1.12.2    2019-04-07 [1] CRAN (R 3.6.1)
 DBI                    1.0.0     2018-05-02 [1] CRAN (R 3.6.1)
 DelayedArray         * 0.10.0    2019-05-02 [1] Bioconductor  
 DelayedMatrixStats     1.6.0     2019-05-02 [1] Bioconductor  
 desc                   1.2.0     2018-05-01 [1] CRAN (R 3.6.1)
 devtools             * 2.1.0     2019-07-06 [1] CRAN (R 3.6.1)
 digest                 0.6.20    2019-07-04 [1] CRAN (R 3.6.1)
 dplyr                * 0.8.3     2019-07-04 [1] CRAN (R 3.6.1)
 DT                   * 0.7       2019-06-11 [1] CRAN (R 3.6.1)
 edgeR                * 3.26.5    2019-06-21 [1] Bioconductor  
 evaluate               0.14      2019-05-28 [1] CRAN (R 3.6.1)
 fitdistrplus           1.0-14    2019-01-23 [1] CRAN (R 3.6.1)
 formatR                1.7       2019-06-11 [1] CRAN (R 3.6.1)
 fs                     1.3.1     2019-05-06 [1] CRAN (R 3.6.1)
 future                 1.14.0    2019-07-02 [1] CRAN (R 3.6.1)
 future.apply           1.3.0     2019-06-18 [1] CRAN (R 3.6.1)
 gbRd                   0.4-11    2012-10-01 [1] CRAN (R 3.6.1)
 gdata                  2.18.0    2017-06-06 [1] CRAN (R 3.6.1)
 GenomeInfoDb         * 1.20.0    2019-05-02 [1] Bioconductor  
 GenomeInfoDbData       1.2.1     2019-07-17 [1] Bioconductor  
 GenomicRanges        * 1.36.0    2019-05-02 [1] Bioconductor  
 ggbeeswarm             0.6.0     2017-08-07 [1] CRAN (R 3.6.1)
 ggplot2              * 3.2.0     2019-06-16 [1] CRAN (R 3.6.1)
 ggrepel                0.8.1     2019-05-07 [1] CRAN (R 3.6.1)
 ggridges               0.5.1     2018-09-27 [1] CRAN (R 3.6.1)
 globals                0.12.4    2018-10-11 [1] CRAN (R 3.6.1)
 glue                   1.3.1     2019-03-12 [1] CRAN (R 3.6.1)
 gplots               * 3.0.1.1   2019-01-27 [1] CRAN (R 3.6.1)
 gridExtra              2.3       2017-09-09 [1] CRAN (R 3.6.1)
 gtable                 0.3.0     2019-03-25 [1] CRAN (R 3.6.1)
 gtools                 3.8.1     2018-06-26 [1] CRAN (R 3.6.1)
 hms                    0.5.0     2019-07-09 [1] CRAN (R 3.6.1)
 htmltools              0.3.6     2017-04-28 [1] CRAN (R 3.6.1)
 htmlwidgets            1.3       2018-09-30 [1] CRAN (R 3.6.1)
 httpuv                 1.5.1     2019-04-05 [1] CRAN (R 3.6.1)
 httr                   1.4.0     2018-12-11 [1] CRAN (R 3.6.1)
 ica                    1.0-2     2018-05-24 [1] CRAN (R 3.6.1)
 igraph                 1.2.4.1   2019-04-22 [1] CRAN (R 3.6.1)
 IRanges              * 2.18.1    2019-05-31 [1] Bioconductor  
 irlba                  2.3.3     2019-02-05 [1] CRAN (R 3.6.1)
 jsonlite               1.6       2018-12-07 [1] CRAN (R 3.6.1)
 KernSmooth             2.23-15   2015-06-29 [1] CRAN (R 3.6.1)
 knitr                * 1.23      2019-05-18 [1] CRAN (R 3.6.1)
 labeling               0.3       2014-08-23 [1] CRAN (R 3.6.1)
 later                  0.8.0     2019-02-11 [1] CRAN (R 3.6.1)
 lattice                0.20-38   2018-11-04 [1] CRAN (R 3.6.1)
 lazyeval               0.2.2     2019-03-15 [1] CRAN (R 3.6.1)
 limma                * 3.40.2    2019-05-17 [1] Bioconductor  
 listenv                0.7.0     2018-01-21 [1] CRAN (R 3.6.1)
 lmtest                 0.9-37    2019-04-30 [1] CRAN (R 3.6.1)
 locfit                 1.5-9.1   2013-04-20 [1] CRAN (R 3.6.1)
 lsei                   1.2-0     2017-10-23 [1] CRAN (R 3.6.1)
 magrittr               1.5       2014-11-22 [1] CRAN (R 3.6.1)
 MASS                   7.3-51.4  2019-03-31 [1] CRAN (R 3.6.1)
 Matrix                 1.2-17    2019-03-22 [1] CRAN (R 3.6.1)
 matrixStats          * 0.54.0    2018-07-23 [1] CRAN (R 3.6.1)
 memoise                1.1.0     2017-04-21 [1] CRAN (R 3.6.1)
 metap                  1.1       2019-02-06 [1] CRAN (R 3.6.1)
 mime                   0.7       2019-06-11 [1] CRAN (R 3.6.1)
 munsell                0.5.0     2018-06-12 [1] CRAN (R 3.6.1)
 nlme                   3.1-140   2019-05-12 [1] CRAN (R 3.6.1)
 npsurv                 0.4-0     2017-10-14 [1] CRAN (R 3.6.1)
 pbapply                1.4-1     2019-07-15 [1] CRAN (R 3.6.1)
 pheatmap             * 1.0.12    2019-01-04 [1] CRAN (R 3.6.1)
 pillar                 1.4.2     2019-06-29 [1] CRAN (R 3.6.1)
 pkgbuild               1.0.4     2019-08-05 [1] CRAN (R 3.6.1)
 pkgconfig              2.0.2     2018-08-16 [1] CRAN (R 3.6.1)
 pkgload                1.0.2     2018-10-29 [1] CRAN (R 3.6.1)
 plotly                 4.9.0     2019-04-10 [1] CRAN (R 3.6.1)
 plyr                   1.8.4     2016-06-08 [1] CRAN (R 3.6.1)
 png                    0.1-7     2013-12-03 [1] CRAN (R 3.6.1)
 prettyunits            1.0.2     2015-07-13 [1] CRAN (R 3.6.1)
 processx               3.4.1     2019-07-18 [1] CRAN (R 3.6.1)
 progress               1.2.2     2019-05-16 [1] CRAN (R 3.6.1)
 promises               1.0.1     2018-04-13 [1] CRAN (R 3.6.1)
 ps                     1.3.0     2018-12-21 [1] CRAN (R 3.6.1)
 purrr                  0.3.2     2019-03-15 [1] CRAN (R 3.6.1)
 R.methodsS3            1.7.1     2016-02-16 [1] CRAN (R 3.6.1)
 R.oo                   1.22.0    2018-04-22 [1] CRAN (R 3.6.1)
 R.utils                2.9.0     2019-06-13 [1] CRAN (R 3.6.1)
 R6                     2.4.0     2019-02-14 [1] CRAN (R 3.6.1)
 RANN                   2.6.1     2019-01-08 [1] CRAN (R 3.6.1)
 RColorBrewer           1.1-2     2014-12-07 [1] CRAN (R 3.6.1)
 Rcpp                   1.0.1     2019-03-17 [1] CRAN (R 3.6.1)
 RCurl                  1.95-4.12 2019-03-04 [1] CRAN (R 3.6.1)
 Rdpack                 0.11-0    2019-04-14 [1] CRAN (R 3.6.1)
 remotes                2.1.0     2019-06-24 [1] CRAN (R 3.6.1)
 reshape2               1.4.3     2017-12-11 [1] CRAN (R 3.6.1)
 reticulate             1.12      2019-04-12 [1] CRAN (R 3.6.1)
 rlang                  0.4.0     2019-06-25 [1] CRAN (R 3.6.1)
 rmarkdown              1.14      2019-07-12 [1] CRAN (R 3.6.1)
 ROCR                   1.0-7     2015-03-26 [1] CRAN (R 3.6.1)
 rprojroot              1.3-2     2018-01-03 [1] CRAN (R 3.6.1)
 RSQLite                2.1.1     2018-05-06 [1] CRAN (R 3.6.1)
 rsvd                   1.0.1     2019-06-02 [1] CRAN (R 3.6.1)
 Rtsne                  0.15      2018-11-10 [1] CRAN (R 3.6.1)
 S4Vectors            * 0.22.0    2019-05-02 [1] Bioconductor  
 scales                 1.0.0     2018-08-09 [1] CRAN (R 3.6.1)
 scater               * 1.12.2    2019-05-24 [1] Bioconductor  
 sctransform            0.2.0     2019-04-12 [1] CRAN (R 3.6.1)
 SDMTools               1.1-221.1 2019-04-18 [1] CRAN (R 3.6.1)
 sessioninfo            1.1.1     2018-11-05 [1] CRAN (R 3.6.1)
 Seurat               * 3.0.2     2019-06-14 [1] CRAN (R 3.6.1)
 shiny                  1.3.2     2019-04-22 [1] CRAN (R 3.6.1)
 SingleCellExperiment * 1.6.0     2019-05-02 [1] Bioconductor  
 stringi                1.4.3     2019-03-12 [1] CRAN (R 3.6.1)
 stringr                1.4.0     2019-02-10 [1] CRAN (R 3.6.1)
 SummarizedExperiment * 1.14.0    2019-05-02 [1] Bioconductor  
 survival               2.44-1.1  2019-04-01 [1] CRAN (R 3.6.1)
 testthat               2.2.1     2019-07-25 [1] CRAN (R 3.6.1)
 tibble                 2.1.3     2019-06-06 [1] CRAN (R 3.6.1)
 tidyr                  0.8.3     2019-03-01 [1] CRAN (R 3.6.1)
 tidyselect             0.2.5     2018-10-11 [1] CRAN (R 3.6.1)
 tsne                   0.1-3     2016-07-15 [1] CRAN (R 3.6.1)
 usethis              * 1.5.1     2019-07-04 [1] CRAN (R 3.6.1)
 vctrs                  0.2.0     2019-07-05 [1] CRAN (R 3.6.1)
 vipor                  0.4.5     2017-03-22 [1] CRAN (R 3.6.1)
 viridis                0.5.1     2018-03-29 [1] CRAN (R 3.6.1)
 viridisLite            0.3.0     2018-02-01 [1] CRAN (R 3.6.1)
 withr                  2.1.2     2018-03-15 [1] CRAN (R 3.6.1)
 xfun                   0.8       2019-06-25 [1] CRAN (R 3.6.1)
 XML                    3.98-1.20 2019-06-06 [1] CRAN (R 3.6.1)
 xtable                 1.8-4     2019-04-21 [1] CRAN (R 3.6.1)
 XVector                0.24.0    2019-05-02 [1] Bioconductor  
 yaml                   2.2.0     2018-07-25 [1] CRAN (R 3.6.1)
 zeallot                0.1.0     2018-01-28 [1] CRAN (R 3.6.1)
 zlibbioc               1.30.0    2019-05-02 [1] Bioconductor  
 zoo                    1.8-6     2019-05-28 [1] CRAN (R 3.6.1)

[1] /home/imallona/soft/R/R-3.6.1/library